% Function : Forward Search Algorithm for Feature Selection
% Parameter:
%   Y: target data that we are trying to predict
%   X: training set, each row is one training example, each column is one feature
%   u: upper limit of selected features, search depth
% Return:
%   I : indexes of selected features
%   R : average residuals of selected features
%   RR: R^2 statistic of selected features
function [I, R, RR] = Forward(Y, X, u)
    if nargin ~= 3
        error('Usage: Forward(Y, X, u)');
    end

    [m, n] = size(X);   % m: number of samples, n: number of features

    % initialize collection
    F  = ones(m, u + 1);    % the training set
    I  = zeros(1, u);       % index of collected features
    R  = zeros(1, u);       % average residuals
    RR = zeros(1, u);       % R^2 statistic

    for k = 1 : u
        r  = zeros(1, n);
        rr = zeros(1, n);

        % search feature i to n
        for i = 1 : n
            % feature i is already in collection?
            if any(i == I)
                continue;
            end

            F(:, k + 1) = X(:, i);
            [~, ~, t, ~, stats] = regress(Y, F(:, 1 : k + 1), 0.05);
            r(i)  = mean(t .^ 2);
            rr(i) = stats(1);
        end

        % select the feature that has the max R^2 statistic
        maxi = find(rr == max(rr));
        maxi = maxi(1);

        % udpate collection
        F(:, k + 1) = X(:, maxi);
        I(k)  = maxi;
        R(k)  = r(maxi);
        RR(k) = rr(maxi);
    end
end
